RNA-Seq: differential expression analysis of controls (PNW8 vs PNW20)

Libraries required

library(limma)
library(edgeR)
library(plgINS)
library(sva)
library(ggplot2)
library(dplyr)
library(patchwork)
library(SummarizedExperiment)
library(reshape2)
library(autoplotly)
library(pheatmap)
library(viridis)
library(RColorBrewer)
library(ggplotify)

Load salmon objects

load("input/SC_controls_rnaseq_salmon.tds.RData")
lab <- salmon

load("input/SC_controls_rnaseq_lit_salmon.tds.RData")
lit <- salmon

counts_lab <- lab@tx.counts
counts_lit <- lit@tx.counts

counts_lab <- counts_lab[,grep(pattern = "Adult", x = colnames(counts_lab))]
counts_lab <- data.frame(Genes = rownames(counts_lab), counts_lab)

counts_lit <- counts_lit[,grep(pattern = "PNW8", x = colnames(counts_lit)), drop = FALSE]
counts_lit <- data.frame(Genes = rownames(counts_lit), counts_lit)

counts <- merge(counts_lab, counts_lit)
rownames(counts) <- counts$Genes

counts <- counts[,-1]

pdata <- data.frame(Samples = colnames(counts), Group = c(rep("PNW20", 6), "PNW8"))

Differeitial analysis using limma

PNW8 vs PNW20

design <- model.matrix(~ 0 + Group, data = pdata)
colnames(design) <- gsub(pattern = "Group", replacement = "", x = colnames(design))

y <- DGEList(counts = counts)
keep <- filterByExpr(y, design, min.count = 10)
y <- y[keep, ]

dds <- calcNormFactors(y)
v <- voom(dds, design = design)

contrast.matrix <- makeContrasts(PNW20 - PNW8, levels = design)
fit <- lmFit(v)

fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

pnw8.pnw20 <- as.data.frame(topTable(fit2, coef = 1, number = Inf))
pnw8.pnw20 <- data.frame(Genes = rownames(pnw8.pnw20), pnw8.pnw20, stringsAsFactors = F)
pnw8.pnw20_sig <- pnw8.pnw20[abs(pnw8.pnw20$logFC) >= 1 & pnw8.pnw20$adj.P.Val <= 0.05, ]

cpm Counts data and pData

y <- DGEList(counts = counts)
dds <- calcNormFactors(y)
cpm <- cpm(dds, log = T)

data <- list(cpm = cpm, pData = pdata)

save(data,
  file = "./output/data_pData.RData", compress = T,
  compression_level = 3
)

Results

dea.list <- list(
  `pnw20 vs pnw8` = as.DEA(pnw8.pnw20)
)

dea.limma <- list(
  `pnw20 vs pnw8` = pnw8.pnw20
)

Save RData files

save(dea.list,
  file = "./output/dea_SC_Controls_pnw8_pnw20.DEA.RData", compress = T,
  compression_level = 3
)

save(dea.limma,
  file = "./output/limma_SC_Controls_pnw8_pnw20.RData", compress = T,
  compression_level = 3
)

writexl::write_xlsx(x = dea.limma, path = "output/dea_results_pnw8_pnw20.xlsx", col_names = T, format_headers = T)

MA plots

res <- pnw8.pnw20
res$threshold <- as.factor(res$adj.P.Val < 0.05)

p1 <- ggplot(data = res, aes(
  x = res$AveExpr,
  y = res$logFC,
  colour = threshold
)) +
  geom_point(alpha = 0.5, size = 1.8) +
  geom_hline(aes(yintercept = 0), colour = "blue", size = 1) +
  ylim(c(
    -ceiling(max(abs(res$logFC))),
    ceiling(max(abs(res$logFC)))
  )) +
  ggtitle("PNW20 vs PNW8") +
  labs(subtitle = "loess fit") +
  xlab("Mean expression") +
  ylab("Log2 Fold Change") +
  theme(
    plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
    plot.subtitle = element_text(size = 16, hjust = 0.5, face = "italic", color = "blue"),
    axis.title.x = element_text(face = "bold", size = 15),
    axis.text.x = element_text(face = "bold", size = 12),
    legend.title = element_text(face = "bold", size = 15),
    legend.text = element_text(size = 14)
  ) +
  scale_colour_discrete(name = "p.adjusted < 0.05") +
  stat_smooth(se = FALSE, method = "loess", color = "red", formula = y ~ x, size = 1)
p1
## Warning: Use of `res$AveExpr` is discouraged. Use `AveExpr` instead.
## Warning: Use of `res$logFC` is discouraged. Use `logFC` instead.
## Warning: Use of `res$AveExpr` is discouraged. Use `AveExpr` instead.
## Warning: Use of `res$logFC` is discouraged. Use `logFC` instead.

Heatmap of most variable genes

xvar <- apply((data$cpm + 1), 1, var)
genes500 <- head(sort(xvar, decreasing = TRUE), n = 500)
x500 <- data$cpm[rownames(data$cpm) %in% names(genes500), ]

pheatmap(mat = x500, scale = "row", show_rownames = F, col = viridis(100), 
         show_colnames = T,
         main = "Clustering of samples based on top 500 variable genes")

PCA

pca <- prcomp(t(data$cpm))
# ggplotly(
  autoplot(pca,
  data = data$pData, colour = "Group",
  frame = TRUE, frame.type = "norm", size = 2
) +
  ggtitle("") + 
    theme(axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),  
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        plot.title = element_text(size = 25, face = "bold"),
        legend.text=element_text(size=20),
        legend.title=element_text(size=20, face = "bold"))
## Too few points to calculate an ellipse

Volcano plots

tab <- lab@tx.info
tab <- tab[,c(4,6,8,12)]
rownames(tab) <- NULL
tab <- tab[!duplicated(tab),]

tab <- merge(tab, pnw8.pnw20, by.x = "transcript", "Genes")

tab$Significant <- "Not significant"
tab$Significant[tab$adj.P.Val <= 0.05 & abs(tab$logFC) >= 1] <- "PNW8"
tab$Significant[tab$Significant == "PNW8"] <- ifelse(tab$logFC[tab$Significant == "PNW8"] > 0, "PNW20", "PNW8")

tab$Significant[tab$Significant == "PNW8"] <- "Ribo0"
tab$Significant[tab$Significant == "PNW20"] <- "polyA"


a <- ggplot(tab, aes(x = logFC, y = -log10(adj.P.Val))) +
  geom_point(aes(color = Significant), size = 2, alpha = 0.5) +
  scale_color_manual(values = c("black", "#0072B5FF", "#BC3C29FF")) +
  xlim(
    -(ceiling(max(abs(tab$logFC)))),
    ceiling(max(abs(tab$logFC)))
  ) +
  ylim(0, max(-log10(tab$adj.P.Val), na.rm = TRUE) + 1) +
  xlab(bquote(~ Log[2] ~ "FC")) +
  ylab(bquote(~ -Log[10] ~ adjusted ~ italic(P))) +
  labs(
    title = "PNW20 vs PNW8"
  ) +
  theme_bw(base_family = "Helvetica") +
  geom_vline(
    xintercept = c(-1, 1),
    linetype = "longdash",
    colour = "black",
    size = 0.4
  ) +
  geom_hline(
    yintercept = -log10(0.05),
    linetype = "longdash",
    colour = "black",
    size = 0.4
  )

a1 <- a + facet_wrap(~type, scales = "free_y")
a2 <- a + facet_wrap(~type2, scales = "free_y")
a

a1

a2

References

report::cite_packages(session = sessionInfo())
## Warning in citation(pkg_name): no date field in DESCRIPTION file of package
## 'plgINS'

SessionInfo

devtools::session_info() %>%
  details::details()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value                       
 version  R version 3.6.2 (2019-12-12)
 os       Ubuntu 16.04.7 LTS          
 system   x86_64, linux-gnu           
 ui       X11                         
 language (EN)                        
 collate  en_US.UTF-8                 
 ctype    en_US.UTF-8                 
 tz       Europe/Zurich               
 date     2020-12-11                  

─ Packages ───────────────────────────────────────────────────────────────────
 package              * version    date       lib
 annotate               1.64.0     2019-10-29 [1]
 AnnotationDbi          1.48.0     2019-10-29 [1]
 assertthat             0.2.1      2019-03-21 [1]
 autoplotly           * 0.1.2      2018-04-21 [1]
 backports              1.1.8      2020-06-17 [1]
 base64enc              0.1-3      2015-07-28 [1]
 Biobase              * 2.46.0     2019-10-29 [1]
 BiocGenerics         * 0.32.0     2019-10-29 [1]
 BiocManager            1.30.10    2019-11-16 [1]
 BiocParallel         * 1.20.1     2019-12-21 [1]
 Biostrings             2.54.0     2019-10-29 [1]
 bit                    4.0.4      2020-08-04 [1]
 bit64                  4.0.2      2020-07-30 [1]
 bitops                 1.0-6      2013-08-17 [1]
 blob                   1.2.1      2020-01-20 [1]
 bookdown               0.20       2020-06-23 [1]
 callr                  3.4.3      2020-03-28 [1]
 checkmate              2.0.0      2020-02-06 [1]
 cli                    2.0.2      2020-02-28 [1]
 cluster                2.1.0      2019-06-19 [1]
 colorspace             1.4-1      2019-03-18 [1]
 crayon                 1.3.4      2017-09-16 [1]
 data.table             1.13.0     2020-07-24 [1]
 DBI                    1.1.0      2019-12-15 [1]
 DelayedArray         * 0.12.3     2020-04-09 [1]
 desc                   1.2.0      2018-05-01 [1]
 DESeq2                 1.26.0     2019-10-29 [1]
 devtools               2.3.1      2020-07-21 [1]
 digest                 0.6.25     2020-02-23 [1]
 dplyr                * 1.0.1      2020-07-31 [1]
 edgeR                * 3.28.1     2020-02-26 [1]
 ellipsis               0.3.1      2020-05-15 [1]
 evaluate               0.14       2019-05-28 [1]
 fansi                  0.4.1      2020-01-08 [1]
 farver                 2.0.3      2020-01-16 [1]
 foreign                0.8-76     2020-03-03 [1]
 Formula                1.2-3      2018-05-03 [1]
 fs                     1.5.0      2020-07-31 [1]
 genefilter           * 1.68.0     2019-10-29 [1]
 geneplotter            1.64.0     2019-10-29 [1]
 generics               0.0.2      2018-11-29 [1]
 GenomeInfoDb         * 1.22.1     2020-03-27 [1]
 GenomeInfoDbData       1.2.2      2019-11-18 [1]
 GenomicRanges        * 1.38.0     2019-10-29 [1]
 GEOquery               2.54.1     2019-11-18 [1]
 ggfortify              0.4.10     2020-07-27 [1]
 ggplot2              * 3.3.2      2020-06-19 [1]
 ggplotify            * 0.0.5      2020-03-12 [1]
 glue                   1.4.1      2020-05-13 [1]
 gridExtra              2.3        2017-09-09 [1]
 gridGraphics           0.5-0      2020-02-25 [1]
 gtable                 0.3.0      2019-03-25 [1]
 Hmisc                  4.4-1      2020-08-10 [1]
 hms                    0.5.3      2020-01-08 [1]
 htmlTable              2.0.1      2020-07-05 [1]
 htmltools              0.5.0      2020-06-16 [1]
 htmlwidgets            1.5.1      2019-10-08 [1]
 httr                   1.4.2      2020-07-20 [1]
 insight                0.9.0      2020-07-20 [1]
 IRanges              * 2.20.2     2020-01-13 [1]
 jpeg                   0.1-8.1    2019-10-24 [1]
 jsonlite               1.7.0      2020-06-25 [1]
 knitr                  1.29       2020-06-23 [1]
 labeling               0.3        2014-08-23 [1]
 lattice                0.20-41    2020-04-02 [1]
 latticeExtra           0.6-29     2019-12-19 [1]
 lazyeval               0.2.2      2019-03-15 [1]
 lifecycle              0.2.0      2020-03-06 [1]
 limma                * 3.42.2     2020-02-03 [1]
 locfit                 1.5-9.4    2020-03-25 [1]
 magrittr               1.5        2014-11-22 [1]
 Matrix                 1.2-18     2019-11-27 [1]
 matrixStats          * 0.56.0     2020-03-13 [1]
 memoise                1.1.0.9000 2020-05-06 [1]
 mgcv                 * 1.8-31     2019-11-09 [1]
 munsell                0.5.0      2018-06-12 [1]
 nlme                 * 3.1-148    2020-05-24 [1]
 nnet                   7.3-14     2020-04-26 [1]
 patchwork            * 1.0.1      2020-06-22 [1]
 pheatmap             * 1.0.12     2019-01-04 [1]
 pillar                 1.4.6      2020-07-10 [1]
 pkgbuild               1.1.0      2020-07-13 [1]
 pkgconfig              2.0.3      2019-09-22 [1]
 pkgload                1.1.0      2020-05-29 [1]
 plgINS               * 0.1.5      2020-07-10 [1]
 plotly                 4.9.2.1    2020-04-04 [1]
 plyr                   1.8.6      2020-03-03 [1]
 png                    0.1-7      2013-12-03 [1]
 prettyunits            1.1.1      2020-01-24 [1]
 processx               3.4.3      2020-07-05 [1]
 ps                     1.3.3      2020-05-08 [1]
 purrr                  0.3.4      2020-04-17 [1]
 R6                     2.4.1      2019-11-12 [1]
 RColorBrewer         * 1.1-2      2014-12-07 [1]
 Rcpp                   1.0.5      2020-07-06 [1]
 RCurl                  1.98-1.2   2020-04-18 [1]
 readr                  1.3.1      2018-12-21 [1]
 remotes                2.2.0      2020-07-21 [1]
 report                 0.1.0      2020-03-19 [1]
 reshape2             * 1.4.4      2020-04-09 [1]
 rlang                  0.4.7      2020-07-09 [1]
 rmarkdown              2.3        2020-06-18 [1]
 rmdformats             0.4.0      2020-06-07 [1]
 rpart                  4.1-15     2019-04-12 [1]
 rprojroot              1.3-2      2018-01-03 [1]
 RSQLite                2.1.4      2019-12-04 [1]
 rstudioapi             0.11       2020-02-07 [1]
 rvcheck                0.1.8      2020-03-01 [1]
 S4Vectors            * 0.24.4     2020-04-09 [1]
 scales                 1.1.1      2020-05-11 [1]
 sessioninfo            1.1.1      2018-11-05 [1]
 SRAdb                  1.48.2     2019-12-24 [1]
 stringi                1.4.6      2020-02-17 [1]
 stringr                1.4.0      2019-02-10 [1]
 SummarizedExperiment * 1.16.1     2019-12-19 [1]
 survival               3.2-3      2020-06-13 [1]
 sva                  * 3.34.0     2019-10-29 [1]
 testthat               2.3.2      2020-03-02 [1]
 tibble                 3.0.3      2020-07-10 [1]
 tidyr                  1.1.0      2020-05-20 [1]
 tidyselect             1.1.0      2020-05-11 [1]
 usethis                1.6.1      2020-04-29 [1]
 vctrs                  0.3.2      2020-07-15 [1]
 viridis              * 0.5.1      2018-03-29 [1]
 viridisLite          * 0.3.0      2018-02-01 [1]
 withr                  2.2.0      2020-04-20 [1]
 writexl                1.3        2020-05-05 [1]
 xfun                   0.16       2020-07-24 [1]
 XML                    3.99-0.3   2020-01-20 [1]
 xml2                   1.3.2      2020-04-23 [1]
 xtable                 1.8-4      2019-04-21 [1]
 XVector                0.26.0     2019-10-29 [1]
 yaml                   2.2.1      2020-02-01 [1]
 zlibbioc               1.32.0     2019-10-29 [1]
 source                            
 Bioconductor                      
 Bioconductor                      
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 Bioconductor                      
 Bioconductor                      
 CRAN (R 3.6.1)                    
 Bioconductor                      
 Bioconductor                      
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 Bioconductor                      
 CRAN (R 3.6.1)                    
 Bioconductor                      
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 Bioconductor                      
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 Bioconductor                      
 Bioconductor                      
 CRAN (R 3.6.1)                    
 Bioconductor                      
 Bioconductor                      
 Bioconductor                      
 Bioconductor                      
 Github (sinhrks/ggfortify@3f4020d)
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 Bioconductor                      
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 Bioconductor                      
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 Github (r-lib/memoise@4aefd9f)    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 local                             
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 Github (easystats/report@dcdd283) 
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 Github (juba/rmdformats@94cd7a3)  
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 Bioconductor                      
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 Bioconductor                      
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 Bioconductor                      
 CRAN (R 3.6.2)                    
 Bioconductor                      
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.1)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.2)                    
 CRAN (R 3.6.1)                    
 Bioconductor                      
 CRAN (R 3.6.2)                    
 Bioconductor                      

[1] /home/ubuntu/R/x86_64-pc-linux-gnu-library/3.6
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library